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Abstract 

The alkali halides have been studied very frequently for their simple structures and 
interesting properties, for example, the high-pressure induced transition^. There are 
several successful models that can be employed. The most famous model of them is 
the Tosi-Fumi interionic potential^ which is an empirical potential derived from the 
experimental data. Whereas, here we developed a new potential model based on the 
Mobius lattice inversion method^'^, which can be derived directly from the cohesive 
energy curve without any experimental data and is more effective than the ab initial 
calculation. With the Mobius interionic potentials, we calculated the structural 
and elastic properties of RbCl crystal. The results are in good agreement with 
experiments. We also studied the high-pressure induced B1-B2 transition of RbCl 
crystal and estimated approximately the transition-point which is about 1.09GPa. 
Further more, we used this potential model to simulate the RbCl melting with 
molecular dynamics. The calculated melt-point is approximately 990K~995K, close 
to the experimental data. 
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1 Introduction 



During the past 60 years the alkali halides have been well studied for their 
thermodynamic, elastic and structural properties. The Tosi-Fumi potential 
modeP is usually employed. Among those available models that can correctly 
describe the alkali halides, the pairwise interionic potentials always include two 
parts: the long-range Coulomb term and the short-range term. The Coulomb 
term is of the form ZiZj/rij, where Vij represents the distance between the ith 
and the jth ions, Zi and Zj represent their effective charges. While there is 
not a definitive form of the short-range term. Tosi-Fumi model uses the forms 
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as 




+ Aijexp{-Bijrij) 



, which represent dipole-dipole, dipole-quadrupole, and repulsion terms respec- 
tively. However, if we have the short-range potential curves, any compatible 
function forms will be fitted to the potential curves with adaptive parameters. 
In pervious models^ the parameterization of potentials is usually obtained by 
empirical fitting which depends on experimental data. These empirical poten- 
tials are only applicable with certainty over the range of interionic distances 
used in the fitting procedure, which may lead to problems if the potential 
is used in a calculation that accesses distances outside this range. Here we 
find the way based on the Mobius lattice inversion method to avoid these 
problems. Using the lattice inversion method we get the short-range poten- 
tial curves solely from the cohesive energy curves which can be calculated 
with quantum mechanics. It is to say that our short-range potential curves 
are derived without any experimental data. This method will be discussed in 
detail later. Another question is how to define the effective charges of cations 
and anions. One possible choice is to assign cations and anions with formal 
charges as Rb-l-1 and Cl-1. But here we use another way to decide the effec- 
tive charges which is depicted in part II. We consider nonintegral charges are 
more reasonable. This paper is organized as follows: First, we introduce our 
Mobius inverse potential model in detail. Further, some structural and elastic 
properties of Bl-RbCl crystal are calculated and compared with experimental 
data. Second, we study the high-pressure- induced transition from Bl struc- 
ture to B2 structure with our model, the transition point is also providSed. 
Finally, we simulate the melting of Bl-RbCl crystal with molecular dynamics 
employing our model and calculate the correct melting temperature with a 
coexist-phase method. Thus we come to the conclusion that our Mobius in- 
verse potential model is successful in describing the RbCl alkali halide and 
this lattice inversion technique should have a wider perspective. 



2 Mobius inverse pairwise interionic potentials 

The main idea of Mobius lattice inversion method can be described simply as 
follows^'®: For a single element crystal, the cohesive energy can be expressed 
as 



where x is the nearest-neighbor distance, ro(n) is the nth neighbor, and (l){x) 
is the pair potential. By a self-multiplicative process from {bQ^n)}, the {b{n)} 
is formed, a multiplicative closed semi-group. This implies that a lot of virtual 
lattice points are involved, but the corresponding virtual coordination number 
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is zero. In the {&(n)}, for any two integers m and n, there is a sole integer k 
such that b{k) = h{m)h{n). Hence, the equation above can be rewritten as 



E{x) = X/ ^{P')4>ip{n)x) 



oo 



(2) 



n=l 



where 




Then the general equation for the pairwise interatomic potential obtained from 
Mobius inversion can be expressed as 



oo 



(j)[x) = 2 ^ I{n)E{h{n)x) (4) 



n=l 



where I(n) has the characteristics of 
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Thus we get the pair potential just from the cohesive energy curve. But we 
encounter some obstacles when we try to get the intcrionic pair potentials from 
the Bl-RbCl cohesive energy curve. First, wc have to know the effective charges 
of cations and anions so that we can remove the Coulomb energy from the total 
energy and remain only the short-range terms. Then we use Mobius lattice 
inversion to obtain the short-range potentials from the so called short-range- 
energy curve. However, there are still three types of short-range potentials as 
Rb-Rb, Rb-Cl, Cl-Cl, and it is impossible to find a lattice constructed only 
with one kind of ion. To solve this problem we build several virtual structures 
of lattice that may be not real for RbCl crystal. These virtual structures 
have some similar traits with the real Bl-strucure. For example, one of them, 
the B3 structure has the same sublatticc of Rb-Rb and Cl-Cl with the Bl 
structure. That means the B3 structure contains the same contributions of 
cation-cation interaction and anion-anion interaction with the Bl structure. 
When the short-range-energy curve of the Bl-structure is subtracted by that of 
the B3-structure, we clearly get the energy that contains only the contribution 
of cation-anion interaction and its relationship with lattice constant. 

The total energies of Bl, B3, Tl, and B2 structures are calculated with 
CASTEP^'S'9 (Cambridge Serial Total Energy Package). The ultra-soft pseu- 
dopotentials for rubidium and chlorine ions are adopted and the GGA-PW 
method is used to cope the exchange-correlation energy. And the k-mesh points 
over Brillouin zone are generated with parameters 4x4x4 for the biggest 
reciprocal space and 1x1x1 for the smallest one by the Monkhorst-Pack- 
scheme^*^ corresponding to lattice constant a. The energy tolerance for SCF 
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convergence is 2 x 106 eV/atom, and the kinetic energy cutoff for plane wave 
basis set is 260 eV. 



To decide the effective charges we process as follows: Since the short-range 
parts are quickly convergent when the lattice constant increases, the total 
contribution of energy with the lattice constant larger than 10.0 A is almost 
completely from the Coulomb part. Then, by using the Madelung constants 
of Bl and B3-structure^^ we can determine the effective charges of ions from 
the energy difference between the total energies of Bl and B3-structure RbCl 
crystals of large lattice constant. 

After having known the effective charges we may calculate the Coulomb en- 
ergies of different structures employing the Ewald summation^^ and remove 
this part of energy from the total cohesive energy. Then we start the lattice 
inversion from the short-range-energy curve. 

For Bl-structure 

E^kia) = E^Xi^) + E^l{a) + E^l{a) + E, (6) 

where ESR means the short-range parts of the total energy, E and 

are the energy contributions of cation-cation, anion-anion and cation- 
anion respectively, Ei refers to the energy of an isolated ion, and a is the 
lattice constant. All energy terms are averaged to each ion. Further we have 



(7) 
(8) 



eI'M = \y. 0+-(^v'(^ + fc-i)^ + (^ + J-i)^ + (j + fc-i)^) (9) 

Where , 4> , and are the short-range interionic potentials of cation- 
cation, anion- anion and cation- anion respectively. 

For B3-structure 

E^l,{a) = E^lia) + E^t{a) + £;fi(a) + E^ (10) 



and 



ElX{a) = \ E <P^A'^^{^ + kY + {i + 3Y + {3 + kY\ (11) 
^--(«) = i . E 4>-il^|{^ + kY-r{^-rJY + {J + kf\ (12) 
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As we can see, the difference is just between the cation-anion interactions. 
So we have 

= E'^lia) - E'^tia) 



(13) 



- 1 , [1^ + " - \V + (' + J - i)^ + + * - > 

(14) 

Note that the isolated-ion energy is also gone. 
This equation may be rewritten as 

1 CXD 

AEl'lix) = -Y: R+-{n)<^+-[B+-{n)x] (15) 

^ n=l 

in which we substitute lattice constant a with the nearest cation-anion distance 
X. Then we follow the way discussed at the beginning of this part to obtain 
the curve of $_)__(a;). The derivation of cation-cation and anion-anion short- 
range potentials is of the same manner. For instance, in order to extract the 
anion-anion short-range interaction we build a Tl-structure lattice. The Tl- 
structure lattice has the same sublattice of Rb-Rb with Bl-structure. So the 
short-range energy difference between those two structures is dedicated by the 
cation-anion and the anion-anion interactions and we can exclude the cation- 
anion part since the cation-anion interionic potential is already known. For the 
cation-cation one we simply use the B2 and Bl-structure lattices, extract the 
cation-cation interaction with the cation-anion and the anion-anion potentials 
available. The whole process can be illustrated in Fig.l. After all the short- 
range-potential curves are attained, we use several forms of functions to fit 
these curves and get the suited parameters. We choose an exponential repulsive 
function for cation-anion and a Morse-stretch function for anion-anion. We find 
that among the three short-range potentials the cation-anion one is excessively 
more important than the other two while the cation-cation one is much smaller 
even than the anion-anion's, almost can be neglected. This leads us to ignore 
the cation-cation short-range interaction for saving calculation time. 
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Rb-Cl 



Cl-Cl 



Effective Cfiarges 



1+ 



q- 



1.8140 



2.4470 6.7525 0.1508 



3.8632 8.5640 0.96930e -0.96930e 



Table 1 

Parameters of potential functions 

So the potentials can be expressed as 
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$ (x) 



D. 



1 — exp 
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(17) 
(18) 



The potential parameters for RbCl are listed in Table I. 



With this potential model we calculate some structural and elastic properties 
of Bl-RbCl crystal at OK. The results are listed in Table II. Compared with 
experimental data at room temperature^^, the results show a good agreement. 



3 High-pressure-induced B1-B2 transition 



Most alkali halides crystallize in the Bl (NaCl-like) structure under ambient 
temperature and pressure but turn into the B2 (CsCl-like) structure under a 
highly external hydrostatic pressure. Since Slater^ first described this phase 
transition in 1924, it has been considered to be one of the simplest first-order 
transitions of alkali halides. However, its mechanism still remains uncertain^^. 
If we restrict this transition to a single-step process without intermediate 
structures'^, there are some mechanisms that may be taken into account. One 
of them is the Buerger mechanism'^, according to which pressure contracts 
the Bl primitive rhombohedral cell along the [111] axis, leading to the B2 
phase. Another one is proposed as the WTM mechanism'^ that concerns with 
the relative displacement of consecutive [100] planes. The third mechanism'^ is 
closely related to the Buerger one and concerns with the orientational relations 
in a pressure-induced transition. 

In our present work, we do not want to discuss the mechanism so we just select 
the Buerger mechanism for its simplicity. We try to explain the reason of tran- 
sition using the energy minimization technique. Here we introduce the Gibbs 
free energy G=U-|-PV, to represent the total crystal energy, where U refers to 
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the whole cohesive energy (inchiding the Coulomb part and the short-range 
parts), P is the external pressure, V is the lattice volume, and PV indicates the 
contribution of external pressure. We don't consider the temperature factor, 
so that all the calculations are under absolute zero. 

Due to its too many freedom-degrees, the Gibbs free energy surface is a hy- 
persurface that can not be drawn on a 2-dimension plane. However, we may 
fix the symmetry of the primitive cell as R3m , change the lattice constant 
simultaneously with the rhombohedral angle along the B1-B2 transition path, 
as described in the Buerger mechanism. Then make the lattice constant as 
X-coordinate, the rhombohedral angle as Y-coordinate and the Gibbs free en- 
ergy as the function of both lattice constant and rhombohedral angle. Thus 
we can draw the Gibbs free energy surface with a 2-dimension expression. 

The Gibbs free energy surface at zero pressure is shown in Fig. 2 (a). Prom 
the scheme we can see at zero pressure the Bl and B2-structure are both 
local minima which mean stable and obviously the Bl-structure has a lower 
energy position than the B2-structure. When we add an external pressure to 
the cell, the whole energy surface will rise and the Bl minimum rises faster 
relatively to the B2 minimum. But the Bl minimum remains lower than the 
B2 minimum until the external pressure is high enough and then a phase 
transition may occur. For this reason, we can expect the Bl-structure to be a 
more stable structure than the B2-structure under a external pressure beneath 
the transition-pressure. 

In order to find the transition-point we calculate the Gibbs free energies of 
relaxed Bl and B2-structure at different external pressures by energy min- 
imization. Both Gibbs free energies increase with the external pressure and 
they have one crosspoint at about l.OQGPa, beneath this value of pressure, 
the Gibbs free energy of Bl-structure is lower than that of B2-structure, while 
above this value, the situation is just the opposite, as shown in Fig. 3. It indi- 
cates that, at an external pressure above l.OQGPa, the B2-structure is more 
stable than Bl, so the transition-point is approximately l.OOGPa. This esti- 
mated transition-point accords well with the experimental data about lOkbar 
in Ref.l3. 

From the Gibbs free energy surface at the transition pressure, shown in Fig. 2 
(b), we can see the Bl and B2 minima are of the same depth which means 
the Bl and B2-structure are of the same stability at this time. However, it is 
difficult for a Bl-structure to turn into a B2-structure spontaneously by energy 
minimization because there is still an energy barrier between the two minima. 
With the further increase of external pressure the Bl minimum will become 
higher and higher. Finally, the Bl-structure will turn to a saddle point while 
the B2-structure remains a minimum and there is clearly a downhill path from 
Bl to B2, as shown in Fig. 2 (c). Then we may achieve the process from B1-B2 
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using the energy minimization method. 



4 Molecular dynamic simulation of RbCl melting 

To perform our model to the thermal properties of Bl-RbCl we employ a MD 
technique. A detailed description of the molecular dynamics method may be 
found elsewhere^^. In our work, the simulations are performed using the MSI 
dynamics engine^". The successful runs rely on the size of system (number 
of particles A^), size of timestep (At), the total running time (total steps of 
run n steps), as well as the parameters of the temperature and pressure control 
methods. After some test runs we find the correct results can be obtained 
with A^ = 512, At = 5 fs, Ugteps = 10000. And we choose Hoover method to 
control the temperature and Andersen method to control the pressure exactly 
at OGPa. The NPT (constant N is the number of particles, T is the temperature 
and P is pressure) ensemble is adopted. 

First we simulate the thermal expansion of Bl-RbCl below the melting point. 
The volume-temperature curve is shown in Fig.4. We also approximately cal- 
culate the linear coefficient of thermal expansion a at different temperatures, 
hsted in Table 111. Compared with the experimental data^^, the calculated 
values are reasonable. 

However, we encounter a discrepance when we try to decide the melt-point 
by heating the perfect lattice until it melts. The abrupt change of volume 
which indicates the solid-liquid transition occius at the temperature of about 
1260K, much higher than the real melting temperature Tm 988K^^. And the 
lattice remains solid at the temperature under which it is expected to have 
been melted. This over-heat problem has been discussed in some papers^^'^^'^^. 
It is said that the calculated melt-point with a perfect lattice as the initial 
configuration will be somewhat larger than the real melt-point, even larger by 
an amount of the order of 20-30%^^. And it is caused by the lack of nucleation 
sites for the liquid phase as the crystal is heated. 

To avoid the over-heat problem and get the correct melt-point we use a coexist- 
phase lattice as the initial configuration^^'^^, which means half of the particles 
are solid and the rest arc liquid. This cocxist-phase lattice can be obtained 
in the following manner: First, we build a superlattice of 512 particles (256 
Rb-I- and 256 CI- respectively). Then, with half of the particles fixed, a MD 
run is carried out at a high enough temperature to make sure the movable 
part can completely melt. Thus the coexist-phase lattice with a common in- 
terface of the solid part and the liquid part is obtained, as shown in Fig. 6 
(a). Using this coexist-phase lattice as the initial configuration, a series of MD 
runs at different temperatures are performed and the final configurations are 
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saved and analyzed. We find the abrupt change of volume happens at about 
995K, very close to the real melting temperature, as shown in Fig.5. From the 
saved trajectory files we calculate the radial distribution function (RDF) and 
the mean-square displacement (MSD) of the final configurations at different 
temperatures. The abrupt changes of RDF and MSD also indicate the occur- 
rence of melting, as shown in Fig. 7, Fig.8. Further more, we pick up two final 
configurations which represent complete-solidification and complete-melting 
as shown in Fig. 6 (b), (c). We can see the final configurations are disparate 
on two sides of the melt-point even though they have the same initial config- 
urations. For comparison, the V-T curve with a perfect lattice as the initial 
configuration is also shown in Fig.5. The difference between the two calcu- 
lated "melt-point" is about 260K. So we prefer the coexist-phase method as 
the correct method to calculate the melt-point of RbCl. 

The MD simulations above are all carried out at GPa. With an external hy- 
drostatic pressure adding to the superlattice, a more serious over-heat problem 
will have to be taken into account. Detailed discussion can be found in Ref.23. 

Actually, the calculated "melt-point" with an initial configuration of prefect 
lattice is not a denotation of melting, but of the mechanical instability of the 
chosen model. One must be sure not to confuse with these two conceptions. 



5 Conclusion 

Based on the Mobius inverse potential model, we have studied the static struc- 
tural and elastic properties of Bl-RbCl crystal. The calculated value are in 
good agreement with experimental data. The high-pressure-induced B1-B2 
transition of RbCl crystal is also studied. The estimated transition-pressure is 
about 1.09GPa, consistent with experimental data^^. Using the MD technique 
we have simulated the melting of RbCl crystal at OGPa and calculated the 
correct melt-point with a coexist-phase method. 

For these apphcations of the Mobius inverse potentials, we consider it as a re- 
liable model although it is derived just from the cohesive energy curve without 
any empirical data employed. 

However, there are still some defects in our model. First, using the Mobius 
inversion method can only obtain the pairwise potential, so that we must add a 
3-body or more potential when it is needed^"^. Second, due to the infinite terms 
of summation as , a quickly convergent E(x) is demanded. That is the reason 
why we need to remove the Coulomb part from the total cohesive energy while 
dealing with the ionic crystals. 
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Despite these flaws, we can see the simphcity and vahdity of this Mobius 
lattice inversion method. This method has successfully been used to obtain 
the interatomic potentials in rare-earth metals and intermetallic compounds^^, 
and now the interionic potentials of alkali halides. 
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Fig. 1. The derivation of Mobius interionic potentials. Three types of virtual struc- 
ture are constructed: B3, Tl, and B2, for the purpose of extracting one from the 
totally three short-range interactions: cation-anion, anion-anion, and cation-cation. 
Of each energy curve, the long-range Coulomb energy has been pre-subtracted from 
the total energy. 
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Fig. 2. The — K Gibbs free energy surfaces of R3m symmetry at (a) zero, (b) Ptr 
and (c) 30GPo, respectively. 
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Fig. 3. The Gibbs free energies of relaxed Bl and B2-structure RbCl crystals increase 
with external pressure. The two curves have a crosspoint at P = 1.09GPa. 
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Fig. 4. The thermal expansion of perfect Bl-RbCl crystal calculated with molecular 
dynamics. Temperature varies from lOOi^ to 900K. The linear coefficient of thermal 
expansion is defined as (1/ L){dL/dT), where L is the lattice constant, and the 
derivative is calculated as a central derivative with a temperature step of 200^. 
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Fig. 5. The volume-temperature curves of the coexist-phase lattice and the per- 
fect lattice. The abrupt change of volume with the one-phase method occurs at a 
temperature about 260K higher than that of the two-phase method. 
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(a) " 

Fig. 6. The initial configuration of a coexist-phase lattice (a) and possible final 
configurations (b)(980i^), (c)(1000i^) after equilibrating the system at different 
temperatures. Dark circles refer to ions that were initially in the solid phase; gray 
circles refer to ions initially in the liquid phase. The system has almost completely 
solidified as (b); and has almost completely melted as (c). 
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Fig. 7. The radial distribution function (RDF) of the final configurations at different 
temperatures. The abrupt change of RDF indicates the occurrence of melting, (a) 
shows the RDF of cation-cation and (b) shows that of cation-anion. 
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Fig. 8. The mean-square displacement (MSD) of all ions as the function of time 
at different temperatures. The change of MSD curves indicates the occurrence of 
melting: at low temperatures the MSD oscillates near its balance, while it increases 
with time at temperatures above melt-point. 
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Table II. Calculated and experimental propeities of Bl-RbCl ciystal, a Conipass98_02 model is 

Fig. 9. Table II. 
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Fig. 10. Table III. 



18 



